SUBROUTINE SUVDFS
! (C) Copyright 1990- ECMWF.
!
! This software is licensed under the terms of the Apache Licence Version 2.0
! which can be obtained at http://www.apache.org/licenses/LICENSE-2.0.
! In applying this licence, ECMWF does not waive the privileges and immunities
! granted to it by virtue of its status as an intergovernmental organisation
! nor does it submit to any jurisdiction.

!**   *SUBROUTINE* *SUVDFS* INITIALIZES COMMON BLOCK *YOEVDFS*

!      A. BELJAARS   E.C.M.W.F.   26/03/90

!      PURPOSE
!      -------

!           *SUBROUTINE *SUVDFS* INITIALIZES THE CONSTANTS NEEDED BY
!      THE EMPIRICAL STABILITY FUNCTIONS AND SETS UP THE TABLE THAT
!      GIVES THE STABILITY PARAMETER ETA (HEIGHT DEVIDED BY
!      *OBUKHOV LENGTH) AS A FUNCTION OF THE GRADIENT *RICHARDSON NUMBER
!      FOR STABLE SITUATIONS. ALSO THE SECOND DERIVATIVES ARE
!      TABULATED FOR LATER USE BY SPLINE INTERPOLATION.
!           *INIVDFS* CHECKS WETHER THE *PHI* AND *PSI* EXPRESSIONS ARE
!      CONSISTENT (BY MEANS OF NUMERICAL DIFFERENTIATION). IF THE
!      FUNCTIONS ARE NOT CONSISTENT, THE ROUTINE ABORTS.

!      INTERFACE
!      ---------

!           *CALL* *SUVDFS* FROM *SUPHEC*

!      METHOD
!      ------

!           *THE ALGEBRAIC EQUATION TO BE SOLVED IS
!      RI=(PHIH/PHIM**2)*ETA, WHERE *PHIH* AND *PHIM* ARE THE GRADIENT
!      STABILITY FUNCTIONS FOR HEAT AND MOMENTUM RESPECTIVELY.
!           *TO SOLVE THE IMPLICIT ALGEBRAIC EQUATION FOR *ETA* AS A
!      FUNCTION OF *RI*, *NEWTON'S METHOD IS EMPLOYED WITH A FIXED
!      NUMBER OF ITERATIONS. THE NECESSARY DERIVATIVES ARE EVALUATED
!      NUMERICALLY.
!           *AFTER COMPLETION OF THE FUNCTION TABLE, THE SECOND
!      DERIVATIVES ARE DERIVED FOR LATER USE BY THE SPLINE
!      INTERPOLATION.

!      EXTERNALS
!      ---------

!           *STATEMENT FUNCTIONS *PHIMS*, *PHIHS*, *PSIMS* AND *PSIHS*

!      REFERENCE
!      ---------

!           *SEE *PRESS ET AL. (1986; NUMERICAL RECIPES - THE ART OF
!      SCIENTIFIC COMPUTING) FOR DETAILS ON THE SPLINE INTERPOLATION.

!      -----------------------------------------------------------------

USE PARKIND1  ,ONLY : JPIM     ,JPRB, JPRD
USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK, JPHOOK

USE YOEVDFS  , ONLY : JPRITBL  ,RITBL    ,ARITBL   ,RCHBA    ,&
 & RCHBB    ,RCHBC    ,RCHBD    ,RCHB23A  ,RCHBBCD  ,&
 & RCHBCD   ,RCHETA   ,RCHETB   ,RCHETC   ,RCHBHDL  ,&
 & RCDHALF  ,RCDHPI2  ,RIMAX    ,DRITBL   ,DRI26  

IMPLICIT NONE
REAL(KIND=JPRB) :: ZU(JPRITBL)

INTEGER(KIND=JPIM) :: ITMAX, JIT, JJP

REAL(KIND=JPRB) :: ZDRV1, ZDRVN, ZEPS, ZETA, ZFXM,&
 & ZFXP, ZP, ZQN, ZRIB, ZUN, ZETAM, ZETAP
REAL(KIND=JPHOOK) :: ZHOOK_HANDLE

!      INSERT STATEMENT FUNCTIONS (PSI-FUNCTIONS)

#include "fcvdfs.h"

!*     1. INITIALIZE CONSTANTS IN COMMON BLOCK
!         ---------- --------- -- ------ -----

!     1.1 CONSTANTS RELATED TO ETA(RI)-TABLE

IF (LHOOK) CALL DR_HOOK('SUVDFS',0,ZHOOK_HANDLE)
RIMAX   = 10._JPRB
DRITBL  = RIMAX/(JPRITBL-1)
DRI26   = DRITBL**2/6._JPRB
RITBL(1)= 0._JPRB

!     1.2 CONSTANTS FOR THE ELLISON TURNER FUNCTIONS (STABLE)

RCHETA=5._JPRB
RCHETB=4._JPRB
RCHETC=8._JPRB
RCHBHDL = 5._JPRB

!     1.21 CONSTANTS FOR HOLTSLAG AND DEBRUIN FUNCTIONS (STABLE PSI)

RCHBA   = 1._JPRB
RCHBB   = 2.0_JPRB/3._JPRB
RCHBC   = 5._JPRB
RCHBD   = 0.35_JPRB
RCHB23A = (2.0_JPRB/3._JPRB)*RCHBA
RCHBBCD = RCHBB*RCHBC/RCHBD
RCHBCD  = RCHBC/RCHBD

!     1.3 CONSTANTS FOR DYER AND HICKS EXPRESSIONS (UNSTABLE)

RCDHALF = 16._JPRB
RCDHPI2 = 2.0_JPRB*ATAN(1.0_JPRB)

!*    3. LOOP OVER TABLE INDEX
!        ---- ---- ----- -----

DO JJP=2, JPRITBL
  ZRIB=(JJP-1)*DRITBL

!     3.1 INITIAL GUESS

  ZETA=RITBL(JJP-1)
  IF (ZRIB  <  0.5_JPRB) THEN
    ITMAX=5
  ELSE
    ITMAX=3
  ENDIF

!     3.2 NEWTON'S ITERATION LOOP WITH DERIVATIVE FROM FINITE DIFFERENCE

  DO JIT=1, ITMAX
    ZEPS=(ZETA+1.0_JPRB)*0.001_JPRB
    ZETAP=ZETA+ZEPS
    ZETAM=ZETA-ZEPS
    ZFXP=ZRIB-ZETAP*PHIHS(ZETAP)/PHIMS(ZETAP)**2
    ZFXM=ZRIB-ZETAM*PHIHS(ZETAM)/PHIMS(ZETAM)**2
    ZETA=ZETA-ZEPS*(ZFXP+ZFXM)/(ZFXP-ZFXM)
  ENDDO

!     3.3 STORE RESULT

  RITBL(JJP)=ZETA

ENDDO

!     4. COMPUTE SECOND DERIVATIVES FROM TABULATED RESULTS
!        ------- ------ ----------- ---- --------- -------

!     4.1 DERIVATIVE AT RI=0. AND RI=RIMAX

ZDRV1=1._JPRB
ZDRVN=(RITBL(JPRITBL)-RITBL(JPRITBL-1))/DRITBL

ARITBL(1)=-0.5_JPRB
ZU(1)=(3._JPRB/DRITBL)*((RITBL(2)-RITBL(1))/DRITBL-ZDRV1)

!     4.2. DECOMPOSITION LOOP OF TRIDIAGONAL MATRIX

DO JJP=2,JPRITBL-1
  ZP=0.5_JPRB*ARITBL(JJP-1)+2.0_JPRB
  ARITBL(JJP)=-0.5_JPRB/ZP
  ZU(JJP)=(6._JPRB*((RITBL(JJP+1)-RITBL(JJP))/DRITBL &
   & -(RITBL(JJP)-RITBL(JJP-1))&
   & /DRITBL)/(2*DRITBL)-0.5_JPRB*ZU(JJP-1))/ZP  
ENDDO

ZQN=0.5_JPRB
ZUN=(3._JPRB/DRITBL)*(ZDRVN-(RITBL(JPRITBL)-RITBL(JPRITBL-1))/DRITBL)
ARITBL(JPRITBL)=(ZUN-ZQN*ZU(JPRITBL-1))/(ZQN*ARITBL(JPRITBL-1)+1.0_JPRB)

!     4.3 BACK SUBSTITUTION

DO JJP=JPRITBL-1,1,-1
  ARITBL(JJP)=ARITBL(JJP)*ARITBL(JJP+1)+ZU(JJP)
ENDDO

IF (LHOOK) CALL DR_HOOK('SUVDFS',1,ZHOOK_HANDLE)
END SUBROUTINE SUVDFS
